#ifndef __FEDCL_H__
#define __FEDCL_H__
//#include <gl.h>
#include<stdio.h>

#include "vector.h"
#include "quaternion.h"
#include "matrix.h"
#include "hybroid.h"
#include "deformation.h"

#define MAX_NUM_HYBROIDS 50
#define MAX_NUM_CONSTRAINTS 100

/* main.c */
void read_Hybroids(const char *filename);
int main(int argc, char** argv);
void dump_window(int wid, char *filename);
void calc_normal (vector *p1, vector *p2, vector *p3, vector *norm);

/* store_data.c */
int store_superquadric_nodes(hybroid *h, char *n);

/* hybroid.c */
hybroid *make_hybroid(FILE *infile);
void allocate_hybroid(hybroid *h);
void init_rect_uv_coords(hybroid *h);
void assemble_hybroid_rectangular_mesh(hybroid *h);
void printh(hybroid *h);

/* user.c */
//void user_reset();
//void user_init() ;
//void user();
//void set_projection();
//void toggle(String s, int b, int f);

/* superquadric.c */
void allocate_sin_cos_arrays(int usize, int vsize);
void compute_sin_cos_arrays(int usize, int vsize);
void compute_superquadric(hybroid *h);
void compute_superquadric_taper_bend(hybroid *h);
void compute_superquadric_node(node *nde, 
			       double r1, double r2, double r3,
			       double e1, double e2);
void compute_superquadric_node_taper_bend(hybroid *h, node *nde, 
			       double r1, double r2, double r3,
			       double e1, double e2);

/* display.c */
//void draw_cube(double size);
//void draw_xyz_axes(double size);
//void draw_overlay_xyz_axes(double size);
//void draw_vector
//  (vector *v, vector *p, double scale, int vector_color);
//void draw_line_between_points
//  (vector *p1, vector *p2, int line_color);
//void display_range_data(range_datapoint *data, int num_data_points);
//void display_range_data_hybroids(hybroid *hybroids[], int numhybroids);
//void display_range_data_array(unsigned char *array, int xsize, int ysize, 
//			      double xscale, double yscale, double zscale);
//void display_bat_spring(vector *bat_position, vector *bat_node_posit);
//void draw_cube_p(vector p, double size);
//void display_range_data_cube(range_datapoint *data, int num_data_points);
//void draw_elements_flat_shading(hybroid *h);

/* element.c */
void compute_element_stiffness_matrix(hybroid *h);
void compute_element_shape_functions_at_integration_points(hybroid *h);

/* deformation.c */
void compute_deformation(hybroid *h);
void compute_element_deformation(hybroid *h);
void compute_superquadric_deformation_e(hybroid *h);
void compute_superquadric_deformation_e_taper_bend(hybroid *h);
void compute_superquadric_deformation_taper_bend(hybroid *h);
void adapt_superquadric_derivatives_taper_bend(hybroid *h, int pos, vector fc, 
                                                         double *ider);
void decrement_superquadric_derivatives_taper_bend(hybroid *h, double *ider);
void compute_superquadric_deformation_derivative_taper_bend(hybroid *h);
void adapt_deformations_at_constraints(constraint *constraints, 
                                             int numconstraints);

/* constraint_forces.c */
void compute_forces_and_rotation_derivatives(hybroid *h);
void compute_derivatives(hybroid *h);
void compute_derivatives1(hybroid *h);
void compute_translation_derivative(hybroid *h);
void compute_rotation_derivative(hybroid *h);
void compute_deformation_derivative(hybroid *h);
void compute_velocity(hybroid *h, int k, vector *v);
void unrotate_fc(hybroid *h, vector fc, vector *rfc);

void compute_constraint_forces(constraint *constraints, int numconstraints);
/*
void compute_constraint_forces(hybroid *h1, hybroid *h2);
*/
void compute_constraint_forces4(hybroid *h1, hybroid *h2);
void compute_constraint_forces3(hybroid *h1, hybroid *h2);
void compute_constraint_forces1(hybroid *h1, hybroid *h2, hybroid *h3);

void compute_other_constraint_value_11(constraint *constri, constraint *constrj
                     , int mpos, int fnum, double t, double tt, double **m);
void compute_other_constraint_value_10(constraint *constri, constraint *constrj
                     , int mpos, int fnum, double t, double tt, double **m);
void compute_other_constraint_value_01(constraint *constri, constraint *constrj
                     , int mpos, int fnum, double t, double tt, double **m);
void compute_other_constraint_value_00(constraint *constri, constraint *constrj
                     , int mpos, double t, double tt, double **m);
void compute_force_influence_on_h1_constraints(constraint *constri, vector fc,
                                                                     int fnum);
void compute_force_influence_on_h2_constraints(constraint *constri, vector fc,
                                                                     int fnum);
void compute_constraint_value(constraint *constri, int fnum, int mpos,
                                    double t, double tt, double **m);
void assign_constraint_forces_to_constraints(constraint *constraints,
					     int numconstraints, double *z);
void remove_constraint_forces(constraint *constraints, 
                                    int numconstraints);
void add_spring_like_forces(constraint *constraints, 
			          int numconstraints);

/* forces.c */
void zeroize_force(hybroid *h);
void unrotate_force(hybroid *h);
void compute_net_applied_force_and_torque(hybroid *h);
void compute_net_damping_force_and_torque(hybroid *h);
void nodes_in_patch
  (hybroid *h, double u, double v, 
   node **n00, node **n01, node **n10, node **n11);
void distribute_force_on_nodes
  (vector *force, double s, double t, double h, double k,
   node *n00, node *n01, node *n10, node *n11);
void compute_range_data_forces(hybroid *h);
void compute_range_data_forces_exhaustive(hybroid *h);
hybroid *nearest_hybroid
  (hybroid **hybroids, int num_hybroids, vector *point, 
   vector *separation, double *distance);
node *nearest_hybroid_node(hybroid *h, vector *point, 
			   vector *separation, double *distance);
void bat_spring_force(vector *bat_position, hybroid *h, double strength);
void display_bat_position(vector *bat_position, double size);
void tabletop_forces
  (hybroid *h, double tabletop_z_position, double strength);
void limit_external_forces
  (hybroid *h, double force_mag_limit);
void image_interp_grad
  (double x, double y, int *image, int xdim, int ydim,
   double *df_dx, double *df_dy);
void compute_image_forces
  (hybroid *h, double z_threshold);
void mouse_spring_force
  (vector *mouse_position, node *nearest_node, double strength);

/* dynamics.c */
void compute_positions(hybroid *h);
void compute_translation(hybroid *h);
void compute_rotation(hybroid *h);
void compute_rigid_motion(hybroid *h);
void compute_motion(hybroid *h);
void compute_rotation_derivatives_and_rotation_matrix(hybroid *h);

/*linear_LU_system_solver.c */
void lubksb(double **a, int n, int *indx, double b[]);
void ludcmp(double **a, int n, int *indx, int *d);
double *dvector(int nl, int nh);
int *ivector(int nl, int nh);
double **dmatrix(int nrl, int nrh, int ncl, int nch);
void free_dvector(double *v, int nl, int nh);
int fullPiv(double **a, int na, int row[], int col[], double eps);
void fullBsv(double **a, int na, int row[], int col[], double x[], double z[], 
                                            double b[], int mode);

/* read-data.c */
range_datapoint *read_ascii_xyz_data(const char *filename, int *num_data_points,
					 double strength_data_force);
int read_ascii_xyz_data_hybroids(const char *filename, hybroid *hybroids[], 
                                  int numhybroids, double strength_data_force);
int read_ascii_xyz_data_hybroids_noise(const char *filename, hybroid *hybroids[], 
                                  int numhybroids, double strength_data_force);
void initialize_datapoints_uv_to_sphere
  (hybroid *h);
void scale_range_data(range_datapoint *data, int num_data_points,
		      double xscale, double yscale, double zscale);
void compute_range_data_center_of_gravity
  (hybroid *h, vector *cog);
unsigned char *read_obvius_image
  (const char *filename, int *xsize, int *ysize);

/* intersection.c */
void bilinear_patch
  (double s, double t, double h, double k,
   vector *p00, vector *p10, vector *p01, vector *p11,
   vector *f);
double squared_distance_to_bilinear_patch
  (vector *p, 
   double s, double t, double h, double k,
   vector *p00, vector *p10, vector *p01, vector *p11);
void grad_squared_distance_to_bilinear_patch
    (vector *p, 
     double s, double t, double h, double k,
     vector *p00, vector *p10, vector *p01, vector *p11,
     double *fu, double *fv);
void compute_nearest_point_on_bilinear_patch
  (vector *p,
   double *u, double *v,
   double h, double k,
   vector *p00, vector *p10, vector *p01, vector *p11,
   double ftol, int *iter, 
   double *fret);
void superquadric
  (double u, double v, 
   double r1, double r2, double r3, double e1, double e2,
   vector *f);
void superquadric_grad_hack
  (double u, double v, 
   double r1, double r2, double r3, double e1, double e2,
   vector *f);
double squared_distance_to_superquadric
  (vector *p, double u, double v,
   double r1, double r2, double r3, double e1, double e2);
void grad_squared_distance_to_superquadric
    (vector *p, 
     double u, double v,
     double r1, double r2, double r3, double e1, double e2,
     double *fu, double *fv);
void compute_nearest_point_on_superquadric
  (vector *p,
   double *u, double *v,
   double r1, double r2, double r3, double e1, double e2,
   double ftol, int *iter, 
   double *fret);
void compute_nearest_point_on_sphere
  (vector *p,
   double *u, double *v,
   double radius,
   double *fret);

#endif __FEDCL_H__